Dynamical evolution of a doubly- quantized vortex imprinted in a Bose-Einstein Condensate 
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The recent experiment by Y. Shin et al. [Phys. Rev. Lett. 93, 160406 (2004)] on the decay of a doubly 
quantized vortex imprinted in 23 Na condensates is analyzed by numerically solving the Gross-Pitaevskii 
equation. Our results, which are in very good quantitative agreement with the experiment, demonstrate 
that the vortex decay is mainly a consequence of dynamical instability. Despite apparent contradictions, 
the local density approach is consistent with the experimental results. The monotonic increase observed 
in the vortex lifetimes is a consequence of the fact that, for large condensates, the measured lifetimes 
incorporate the time it takes for the initial perturbation to reach the central slice. When considered locally, 
the splitting occurs approximately at the same time in every condensate, regardless of its size. 

PACS numbers: 03.75.Kk, 03.75.Lm, 67.90.+Z 
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Since the creation of the first dilute-gas Bose-Einstein 
condensates there has been great interest in characterizing 
their superfiuid properties. This has stimulated a great deal 
of theoretical and experimental work aimed at studying the 
rotational properties of dilute Bose gases [ 1 ] and, in par- 
ticular, the nucleation and stability properties of vortices 
J2]. Numerous experiments have succeeded in generating 
vortices [3]. In practically all of them, the vorticity appears 
concentrated in a number of singly quantized vortices. This 
is a consequence of the fact that multiply quantized vortices 
are energetically unstable against their splitting in an array 
of vortices with unit topological charge [4]. Multiply quan- 
tized vortices are also dynamically unstable, which implies 
that they decay even in the zero-temperature limit [5]. 

Recently, Leanhardt et al. [6] obtained multiply quan- 
tized vortices in a gaseous Bose-Einstein condensate by 
using a topological phase-imprinting technique proposed 
by Nakahara et al. Q. They started from nonrotating 23 Na 
condensates in the |1, —1 > and |2, +2 > hyperfine states. 
By adiabatically inverting the magnetic bias field, the ini- 
tial states were transformed into vortex states with axial 
angular momentum per particle 2h and —Ah, respectively. 
This has opened the possibility for studying experimentally 
the stability of multiply quantized vortices, and has stimu- 
lated new theoretical works. In particular, Mottonen et al. 
10] have studied numerically the stability of a doubly quan- 
tized vortex in a cylindrical condensate as a function of the 
(dimensionless) interaction strength per unit length along 
the condensate axis, an z , where a is the s-wave scattering 
length and n z = J j ^ (r ) | 2 dx dy is the linear atom density 
along z. They found a series of quasiperiodic instability re- 
gions in the parameter space of an z ]5|]. The first two of 
them (the only ones relevant to this work) correspond to 
an z < 3 and 11.4 < an z < 16, respectively. By com- 
paring with the solution of the Gross-Pitaevskii equation 
for a harmonically trapped cigar-shaped condensate these 
authors conclude that a doubly quantized vortex is dynam- 
ically unstable, and it is in the above instability regions 
where the vortex decay is initiated, giving rise in the pro- 
cess to a pair of intertwining singly-quantized vortices. 



Shortly after, a second experiment was carried out at 
MIT HI, aimed at studying the characteristic time scale of 
the splitting process as a function of the interaction strength 
an z=0 . Experimental results were obtained from integrated 
absorption images along the condensate axis after 15 ms of 
ballistic expansion. A key ingredient was the fact that in or- 
der to increase the visibility of the vortex cores, absorption 
images were restricted to a 30 fim thick central slice of the 
condensate. This experiment showed that, even at T ~ 0, 
doubly quantized vortices decay into two singly quantized 
vortices on a time scale that is longer at higher atom den- 
sity, with an observed lifetime that increases monotonically 
with an z=0 , showing no quasiperiodic behavior. 

These results, pose a number of open questions. First, 
the observed monotonic behavior seems to be not consis- 
tent with the local density approach proposed in Ref J^]. 
According to this proposal, as the atom density increases, 
for an z=0 ~ 12, the second instability region mentioned 
above arises within the central slice. One then would ex- 
pect this fact to manifest as a decrease in the observed life- 
time, which does not occur. This raises some doubts on the 
validity of the local density approach [9]. Second, numer- 
ical work has appeared recently suggesting that the vortex 
decay is a consequence of thermal fluctuations [ 10]. Even 
though it seems reasonable to attribute the splitting to a dy- 
namical instability, it would be desirable to have theoretical 
results in good quantitative agreement with the experiment, 
which could corroborate this assumption. Finally, the de- 
tailed dynamical behavior of the vortex along the entire z- 
axis can be especially relevant for characterizing the split- 
ting process in very elongated condensates such as those 
studied in the MIT experiment. Different regions along the 
z-axis could exhibit different local behavior which could 
be determinant to understand the process. The experimen- 
tal results cannot provide this kind of information, and thus 
it would be desirable to complement them with the detailed 
dynamics of the vortex along the entire condensate. 

In this Letter we address the above questions by perform- 
ing a realistic computer simulation of the MIT experiment. 
Our results, which are in very good quantitative agreement 
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with the experiment, enable us to understand the experi- 
mental results by providing a detailed visualization of the 
splitting process. They confirm that the decay of the vortex 
core is essentially a dissipationless process driven by a dy- 
namical instability. They also indicate that the local density 
picture proposed in Ref. [8] provides the key ingredients to 
interpret properly the splitting process. When considered 
locally, the decay is characterized by the sole parameter 
an z and it is always initiated in those regions along the z- 
axis where an z ~ 1.5 or an z ~ 13.75, on a time scale that 
is in all cases of the order of 15 ms. For large condensates 
(an z=0 ^ 3), the lifetime observed in the central slice in- 
corporates the time it takes for the nonlinear perturbation 
to propagate from the location where the splitting is ini- 
tiated to the final position where it is eventually detected. 
This nonlocal process explains the observed monotonic in- 
crease in the vortex lifetime and makes the local density 
picture compatible with the experimental results. 

In the zero-temperature limit, the dynamics of dilute 
Bose gases is accurately described by the Gross-Pitaevskii 
equation (GPE), which governs the time evolution of the 
condensate wave function 

ih^={-^V 2 + V(v) + gN\^(v,t)\ 2 ^ (1) 

where N is the number of atoms, g = 4:irh 2 a/m is the 
interaction strength, and ^(r) = ^m(uj 2 + col) is the ex- 
ternal confining potential. 

For the 23 Na condensates used in the MIT experiment, 
a = 2.75 nm and u; r /2ir = 220 Hz. The experimen- 
tal results were obtained using, without distinction, three 
different axial trap frequencies, lo z /2tt = 2.7, 3.7, and 
12.1 Hz. Integration of the 3D GPE for such highly elon- 
gated condensates is a very demanding computational task. 
As the system evolves in time it develops a complex fine 
structure that can only be properly resolved by using very 
large basis or gridpoint sets. To verify the convergence and 
accuracy of our numerical results we have implemented 
two different integration methods: a Laguerre-Fourier and 
a Laguerre-Hermite-Fourier pseudospectral method. The 
evolution in time has been carried out by a third-order 
Adams-Bashforth time-marching scheme. The same re- 
sults have been obtained with the two integration methods 
and using very different basis sets and time steps. 

We have solved the GPE starting, in all cases, from the 
stationary state compatible with a doubly quantized vortex. 
A quantum system is dynamically unstable when it is un- 
stable against arbitrarily small perturbations. To determine 
whether the vortex decay can be unambiguously attributed 
to a dynamical instability, at t = 0, we introduce a small 
fluctuation in the trapping potential. Specifically, we intro- 
duce a 1% quadrupolar perturbation for a short interval of 
0.3 ms. Such a small perturbation produces an almost un- 
detectable change in both the energy and the angular mo- 
mentum per particle (AE/E, AL Z /L Z < 10~ 5 ). In order 




FIG. 1 : (color online) (a) In-trap absorption images of a conden- 
sate with an z= Q = 1.48. (b) Same as (a) after 15 ms of ballistic 
expansion. Lengths are in units of the axial trap, a z = 6.05 /im. 
(c) Predicted splitting times as a function of an z= o for three dif- 
ferent perturbations (curves 1-3). Curve 4 shows the excitation 
spectrum of the unstable modes. 



to compare with the experiment we have followed the same 
procedure as used at MIT. Condensates with different val- 
ues of an z=0 were produced by varying N, and integrated 
absorption images of a 30 p thick central slice were then 
generated at different times. In all cases the initial vortex 
decayed into a pair of singly quantized vortices. We have 
also considered the potential effect of dissipative processes 
by introducing a phenomenological imaginary time and, 
for values consistent with the experiment, we have found it 
to be negligible. The vortex lifetime was inferred from the 
absorption images by identifying the instant at which two 
vortex cores can be resolved unambiguously. An example 
is shown in Fig. We have also analyzed the effect of the 
ballistic expansion by solving numerically the correspond- 
ing GPE in a number of representative cases. This process 
modifies the predicted lifetime by only 1 — 2 ms approx- 
imately (Fig. which is a consequence of the fact that 
most of the expansion is a mere scale transformation. Since 
this correction is of the order of measurement uncertainties 
we neglect it in what follows. 

To investigate the effects of both the strength and the 
duration of the perturbation we have considered two ad- 
ditional cases: a 4% and a 1% quadrupolar perturbation 
acting for 0.3 ms and 2 ms, respectively. The correspond- 
ing numerical results, obtained for an axial trap frequency 
v z = 12 Hz, are shown in Fig. Even though the local 
value of an z essentially determines where and when the 
splitting starts, the sole parameter an z=0 is not sufficient 
to characterize completely the observed splitting times in 
the central slice. For instance, for the 4% quadrupolar per- 
turbation the predicted splitting time for a condensate with 
an z= o = 3.27 in the trap with v z = 3.7 Hz turns out to 
be 55 ms, whereas for a condensate with the same value of 
cm z=0 but in the 12 Hz trap one obtains 33 ms. In general, 
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FIG. 2: (color online) (a) Predicted splitting times as a function 
of an z=0 for a 4% quadrupolar perturbation acting during 0.3 
ms. (b) Density isosurface of a condensate with an z= o — 13.75 
at t = 75 ms. The corresponding axial absorption image of the 
central slice (rectangle in the figure) is also shown. Lengths are 
in units of a z — 6.05 fim. 



for an z=0 ^ 3, the more elongated the condensate, the 
longer the measured lifetime. This is a consequence of the 
fact that measured lifetimes incorporate the time it takes for 
the initial perturbation to reach the central slice. Since the 
shortest lifetimes occur for the largest axial frequency (12 
Hz), this is the only case one has to consider: experimental 
data corresponding to two visible cores must lie above the 
theoretical curve. Note that this does not prevent data cor- 
responding to a single core (those obtained with v z = 2.7 
or 3.7 Hz) from also lying above the theoretical curve. 

While the perturbation strength has an important influ- 
ence on the predicted lifetimes, its duration seems to be 
of little importance (Fig. [Q;). The main effect of stronger 
perturbations is to shift the predicted curve toward shorter 
times. The best agreement with the experiment is found for 
the 4% quadrupolar perturbation (curve 1). However, the 
predicted lifetimes are somewhat longer than the experi- 
mental ones. Note that the experimental results do not in- 
clude the 12 ms spent in the inversion of the axial magnetic 
field B z (the vortex imprinting). Since the vortex already 
forms as B z = 111 ill , from this process (the preparation 
time) one reasonably can expect a contribution of about 

5 — 6 ms. Fig. |2]shows the corresponding theoretical curve 
(from which we have subtracted a conservative amount of 

6 ms to account for the preparation time) along with the 
experimental data of the MIT experiment. The good agree- 
ment demonstrates that the decay is mainly driven by a dy- 
namical instability. 

During the vortex preparation, as B z vanishes, all of the 
three F = 1 components appear in the trap. Thus, for a 
short interval around the instant of preparation the spino- 
rial character of the multicomponent BEC cannot be ne- 
glected iflTll . The weak- field seeking state becomes neces- 
sarily perturbed by the other components. In the presence 



of gravitational interaction (perpendicular to the z-axis in 
the experimental setup) each component is shifted from the 
z-axis by a different amount Jl2ll . This introduces an effec- 
tive non-axially symmetric perturbation on the dynamical 
evolution of the weak-field seeking component via the cor- 
responding interaction terms. It seems reasonable to think 
that the quadrupolar component of this complex effective 
perturbation might be responsible for the vortex decay. 

The theoretical curve in Fig. |3 displays a clear mini- 
mum about an z =o — 1.5. It is instructive to put this min- 
imum in relation to the excitation spectrum of the unsta- 
ble (imaginary frequency) modes (see curve 4 in Fig. 

For small condensates with an z=0 < 1.5 all of the 
z-axis lies within the first instability region. For such con- 
densates only a limited subset of the unstable modes asso- 
ciated with that instability region become accessible. The 
smaller the condensate the smaller the maximum unstable 
frequency and, consequently, the longer the observed decay 
time. For condensates with an z=Q > 1.5 all of the unsta- 
ble frequencies become accessible. In such condensates, 
as the number of particles increases the first instability re- 
gion (corresponding to those values along the z-axis where 
an z < 3) moves progressively away from the central slice, 
towards the edges of the condensate, occupying symmetric 
positions around z = 0. As a consequence, the splitting of 
the vortex core has to propagate from those regions where 
it is initiated to the central slice before it can be detected. 
This nonlocal process is responsible for the monotonic in- 
crease in the lifetimes for an z=0 ^ 3, and explains the 
minimum predicted about an z=Q ~ 1.5. According to a 
local density picture, in principle, one also would expect a 
second minimum about an z=0 ~ 13.75. Even though no 
such a minimum occurs, a clear change in the slope of the 
predicted curve can be identify as an z=Q enters the second 
instability region, indicating a different physical behavior. 

In Fig. |3]we have calculated the density along the z-axis 
of every condensate (characterized by the different an z=0 ), 
at the same instant of time (t = 15 ms). This axial den- 
sity has been renormalized in each z-plane in such a way 
that the maximum density in that plane is 1. Thus, for a 
given an z=0 , dark zones in the plot density indicate those 
regions along the condensate axis where the vortex split- 
ting has already begun. As is evident from Fig. |3] when 
considered as a local process, the splitting takes place ap- 
proximately at the same time (on a timescale of the order 
of 15 ms) in every condensate, regardless of its size. In all 
cases the process starts precisely at the predicted instabil- 
ity regions. Fig. |3]also shows that for small condensates 
with an z=0 < 3 the splitting occurs almost simultaneously 
along the entire z-axis, which is a consequence of the fact 
that in such condensates the entire z-axis lies within the 
first instability region. As a result, the initial vortex decays 
into a pair of nearly straight unit-charge vortices. Note that, 
for an z=0 < 1 • 5 the splitting takes somewhat longer times, 
as expected. As the condensate size increases, the vortex 
splitting becomes a localized phenomenon that propagates 
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FIG. 3: (color online) Density along the z-axis as a function of 
an z= o at t = 15 ms. Lengths are in units of a z = 6.05 /im. 



along the z-axis, producing in this case a pair of intertwin- 
ing singly-quantized vortices |8], which is a consequence 
of the dephase in time between different z -planes. For 
an z=0 > 12 a new instability region appears at the center 
of the condensate and the vortex splitting is firstly initiated 
at the edges of the condensate and a few milliseconds later 
at the central slice (a consequence of the smaller value of 
the corresponding maximum unstable frequency). 

Fig. |4] shows the evolution in time of the splitting pro- 
cess for a condensate with an z=0 = 13.75. The first and 
second instability regions correspond to the shaded zones 
in Fig. |4^. This figure clearly shows that, at t = 25 ms, 
the splitting process has already begun in both the edges 
and the center of the condensate. However, the two vortex 
cores overlap (Fig. @J>) and thus cannot be experimentally 
resolved until much longer times. At t ~ 70 ms (Fig. |4|:) 
the vortex cores begin to disentangle in such a way that they 
can be unambiguously resolved at t = 75 ms (Fig. |3p). 
Thus, despite appearances no contradiction exists with the 
local density picture. No decrease is observed in the de- 
cay times due to the fact that the time required to resolve 
the two vortex cores is much longer than that it takes for 
the instability originating at the edges of the condensate to 
reach the central slice. 

In conclusion, we have analyzed the MIT experiment by 
solving the full GPE. Our results, which are in very good 
quantitative agreement with the experiment, demonstrate 
that the vortex decay is mainly driven by a dynamical in- 
stability, and allow a better understanding of the splitting 
process. Despite apparent contradictions, the local den- 
sity approach is consistent with the experimental results. 
The monotonic increase observed in the vortex lifetimes is 
a consequence of the fact that, for large condensates, the 
measured lifetimes incorporate the time it takes for the ini- 
tial perturbation to reach the central slice. When consid- 
ered locally, the splitting occurs approximately at the same 
time in every condensate, regardless of its size. 

This work has been supported by Ministerio de Edu- 
cacion y Ciencia (Spain) and FEDER fund (EU) (contract 
No. Fis2005-02886). 

Note added-After this work was finished, we learned 



(a) 25 ms (b) 40 rns |«1 70 ms 




FIG. 4: (color online) Evolution in time of the splitting process 
of a condensate with an z —o = 13.75. The shaded zones in (a) 
indicate the instability regions. The corresponding axial absorp- 
tion images of the central slice (rectangle in the figure) are also 
shown. Lengths are in units of a z = 6.05 /im. 



about a very recent preprint |quant-ph/0605125| in which 
the same experiment is analyzed. 
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